1. Field
The invention is in the field of the automated analysis and correction of n-dimensional spectral data, such as spectral or imaging data obtained through nuclear magnetic resonance, electron spin resonance, ion cyclotron mass spectrometry, infrared spectrometry, etc.; where n is greater than or equal to 1.
2. State of the Art
A Nuclear Magnetic Resonance (NMR) Spectrometer subjects a sample to a magnetic field, stimulates the sample with an electromagnetic pulse, and measures the responses of the sample. Such a device comprises a magnet to generate a magnetic field, a source of pulsed electromagnetic energy, and a sensor for monitoring the response of a sample subjected to both the magnetic field and to the pulsed electromagnetic energy. An introduction to the principles of NMR spectrometry appears in Chapter 8 of Alan G. Marshall and Francis R. Verdun, Fourier Transforms in NMR, Optical, and Mass Spectrometry, A User's Handbook, Elsevier, 1990. An Electron Spin Resonance (ESR) spectrometer is structurally similar.
An Ion Cyclotron Resonance (ICR) Mass Spectrometer also comprises a magnet for generating a magnetic field, a source of electromagnetic energy, and a sensor for monitoring the response of a sample to the electromagnetic energy. The sample for ICR consists of charged particles, usually molecules, in a vacuum. An introduction to ICR appears in Chapter 7 of Marshall and Verdun (id.)
Infrared and visible spectrometers subject a sample to a beam of electromagnetic energy, and have a detector to measure the energy transmitted through or modified by the sample. The frequency response of the sample may be obtained by subjecting the sample to a scanned frequency source, or by using a broadband source with a scanning interferometer. The essentials of modern infrared and visible transmission spectrometers are described in Chapter 9 of Marshall and Verdun (id.).
Various methods are used to obtain spectral or imaging data for a variety of purposes. The data obtained may be collected in the time or frequency domain. Collected data can be converted between the time and frequency domain by sampling and digitizing a time domain signal and using a digital computing device to perform a Fourier Transform (FT). Data are usually obtained by perturbing an object to be studied and measuring the response from the object. For example, in nuclear magnetic resonance (NMR), an object is subjected to magnetic fields varying in both time and space. In the case of magnetic resonance imaging (MRI), the spatial variation of the fields is of crucial importance in spatially selecting portions of the sample for measurement. In the case of both imaging and spectroscopy, pulses of radio or microwave frequency (RF) electromagnetic energy are applied to cause a resonance phenomenon. Detectors measure the magnetic fluctuation in the object under study and provide data that are evaluated to determine various characteristics of the object. These magnetic fluctuations are oscillatory, and the frequency, phase, relaxation time and amplitude of the oscillations carry information about the sample. Because of various characteristics of the instrumentation, signal properties such as frequencies, phases, lineshapes, and amplitudes may be distorted and in addition may change over extended periods of time. In order to obtain maximum information from the data, the data must either be corrected to remove these instrumental artifacts or be automatically analyzed in a way insensitive to these distortions.
Many spectrometers do not directly digitize the signal obtained by detecting the response of the sample because of the limited bandwidth of available Analog to Digital Converters (ADC's). These spectrometers mix the detected signal with the output of a local oscillator to heterodyne the signal to a lower frequency range in the same manner that a conventional superheterodyne radio receiver shifts signals into its Intermediate Frequency range. In NMR, the local oscillator is known as the "transmitter", and the frequency at which it operates is known as the "transmitter offset."
To stabilize the position of signals, various hardware solutions to lock the field/frequency ratio are in routine use as described in Derek Shaw Fourier Transform N. M. R. Spectroscopy (second edition); Elsevier: Amsterdam, 1984; p. 142ff. Most spectrometers use heteronuclear locking on a lock substance added to the sample, typically keeping the position of one deuterium resonance constant while observing other nuclei such as .sup.1 H or .sup.13 C. Homonuclear locking is used in some continuous wave (c.w.) proton spectrometers. With homonuclear locking, frequently a small amount of a substance is added to the sample for locking purposes, such as tetramethylsiloxane (TMS), containing the same nucleus as the nucleus being observed, but having a local molecular environment such that it produces an easily recognizable resonance. The prior art uses the frequency value corresponding to the experimental signal from the additive to determine the needed correction. Spectra may also be taken without feedback control of the instrument, but the resulting spectra must generally be interpreted with reference to the position of a similar reference signal. Further, if spectra are to be averaged for noise reduction, the spectra must be corrected for frequency drift of the instrument, known as positional shifts.
Addition of an additive chosen to provide a recognizable locking signal is undesirable because this technique requires that the sample be polluted with the additive. The additive might modify the sample characteristics and there may be an overlap between the signal obtained from the additive and the desired signal from the sample being tested.
Sample saturation results from applying on the average more energy to the sample than the sample can dissipate. The generally used method to detect signal saturation is to increase the excitation power and to see whether or not the sample response increases as expected for an unsaturated sample. The inventor is not aware of prior art for the automatic determination of sample saturation through monitoring determined signal parameters nor of an automated way to control sample saturation.
The most common method currently in use to phase correct spectral data is to manually adjust coefficients of a phase correction function until the experimental data exhibit expected visual properties as observed by a spectroscopist. The extent of correction that can be accomplished manually is limited, and the outcome is subject to the biases of the spectroscopist. With most steps of the data collection and processing being automated, manual phasing is a weak link in the process.
Where frequency domain data are used, such data are generally obtained from time domain data through Fourier transformation. Modern spectrometers and imaging devices use quadrature detection and a complex Fourier transform of the acquired data yields a complex spectrum, i.e., a spectrum consisting of a real and an imaginary part. Ideally, the real part of such a spectrum would contain pure absorption mode responses and the imaginary part would contain pure dispersion mode responses. However, a mismatch between reference and detector reference phase introduces frequency-independent phase shifts, while delays between excitation and acquisition, frequency sweep excitation, and delays introduced by electronic filtering produce a frequency dependent mixture of both modes. For higher-dimensional spectra there is a similar problem associated with each dimension, thus manual correction of multidimensional spectra can become extremely tedious or even impossible.
Methods to avoid phase correction by processing and displaying the information in a phase independent way are in frequent use. However, such methods reduce the information content of the data and for most purposes phase sensitive data are preferable.
A pure absorption mode spectrum can be obtained from phase sensitive data by properly mixing real and imaginary parts at every point in the spectrum, a process called phase correction. The pure absorption mode yields data with the best resolution and sensitivity. Also aliased peaks can be identified by their anomalous phases, a detection which is not possible for phase insensitive spectra. In addition, positive and negative intensities can be distinguished and, if the spectra are taken carefully, the area under an absorption mode signal is proportional to the number of nuclei generating the signal.
A properly phased complex one-dimensional (1D) spectrum can be written as EQU F(.omega.)=A(.omega.)+i D(.omega.) [1 ]
wherein the real portion of the spectrum, Re[F(.omega.)]=A(.omega.), contains the absorption mode signals, and the imaginary portion of the spectrum, Im[F(107)]=D(107 ), the dispersion mode signals. In such equations, the angular frequency, .omega., has units of radians and is related to the frequency, by =.omega./(2.pi.). An experimental spectrum, G(.omega.) can show frequency-independent and frequency-dependent phase shifts and is related to the well-phased spectrum F(.omega.) by the phase function, .phi.(.omega.), as follows: EQU G(.omega.)=F(.omega.) e.sup.-.phi.(.omega.) [ 2]
Another notation in frequent use employs the opposite sign in the exponent and thereby shifts the phase angles by .pi.. The relative phase shift between real and imaginary spectra would be -.pi./2 instead of .pi./2. Phase correction attempts to estimate F(.omega.) from the observed G(.omega.), typically by determining an approximate phase function .phi.'(.omega.). Often, for two or more test signals in the spectrum G(.omega.), coefficients of the phase connection function are adjusted manually so that the corrected phase angle at the transition frequency of each of the chosen lines looks optimally adjusted. Underdigitized spectra always appear to be somewhat poorly phased due to the sinc character of the observed lineshape, making it desirable to use more objective criteria for phase adjustment.
Since signal phases are generally a slowly varying continuous function of the resonance frequency, a best-fit polynomial can be used to interpolate the phases of the test lines so as to correct every point in the spectrum. Phases in 1D spectra can only be determined within a multiple of 2.pi.. In higher-dimensional spectra a simultaneous phase shift in two orthogonal directions of .pi. does not change the appearance of the spectrum and in either case, optimum phases for the test resonances might need to be shifted by 2.pi. or .pi. before interpolation. The first two terms of a power series, the zero.sup.th and first order corrections .phi..sub.0, and .phi..sub.1, .phi.respectively, in: EQU .phi.'(.omega.)=.phi..sub.0 +.phi..sub.1 .omega.+.phi..sub.2 .omega..sup.2 +. [3]
are often sufficient to achieve a satisfactory phase adjustment. This approach of "frequency-dependent" phase shifts is the only widely used method but it is not an accurate way of phase correcting a spectrum. In reality, spectra are a sum of lines, each with its own frequency-independent phase. The phase angle of a signal tail is determined not by its position in the spectrum, but by the phase angle at the signal's transition frequency. The use of the approximate low-order phase function for phase correction of the spectrum G(.omega.) according t o EQU F(.omega.)=G(.omega.) e.sup.i .phi.'(.omega.) [4 ]
modifies the lineshape of the signals by the frequency-dependent phase shift of this equation and introduces baseline undulations.
Typically, "difficult" 1D spectra and some 2D spectra are phased manually using the experience of a trained spectroscopist. In higher-dimensional spectra and in imaging applications, "phase-independent" methods are currently often used because the amount of data involved are too large for manual phasing.
Several methods have been proposed for the automated phase correction of 1D spectra. Ernst published the first article on autophasing methods in the Journal of Magnetic Resonance 1969, 1, 7. The first method described by Ernst calculates the zero.sup.th order phase angle .phi..sub.o with the following relationship: ##EQU1## where I.sub.q is the integral of the imaginary part of the entire spectrum and I.sub.r, is the corresponding integral of the real spectrum. The calculated angle, .phi..sub.o, is then used to phase correct the data. By replacing the integrals by amplitudes in Eq. 5, Montigny et al., in a method described in Analytical Chemistry 1990, 62, 864, determine the phases for individual points in the spectrum near the signal maximum. These determined phases are then used to correct the data. The second method described by Ernst is based on the ratio of maximum to minimum signal excursion. This ratio becomes infinite for pure absorption mode signals and it is unity for a dispersion mode signal. Ernst applied phase corrections to the spectrum to maximize this ratio by an iterative procedure. When the ratio is maximized, the data are supposedly phased. Using a Simplex algorithm, Siegel, in a method described in Analytica Chimica Acta, 1981, 133, 103, used the maximization of signal height and the minimization of the area below the baseline as criteria for optimization. Daubenfeld et al., as described in the Journal of Magnetic Resonance 1985, 62, 195, looked at both Ernst's second method criteria and Siegel's criteria and used an interpolation between spectral points with a Lorentzian lineshape model to improve the optimization based on the criteria of maximum peak area, maximum signal height, and a new criterion, minimum remaining phase deviation, i.e., phase angle. If any one of these criteria are met, the data are assumed to be correctly phased by Daubenfeld's method.
Brown et al. Journal of Magnetic Resonance 1989, 85, 15, used the criteria of flat baseline and narrow linewidth at the base of the signal for phasing. Maximization of the number of spectral points inside a small region of amplitudes defined by the noise level around a flat baseline will yield phased spectra with signals in positive or negative absorption. Van Vaals and van Gerwen Journal of Magnetic Resonance 1990, 86, 127, proposed to determine the best spectral phasing by iteratively recalculating a crude spectral model with varying phase distortions, determining the phase function from the model and fitting this function to the phase function of the measured spectrum over both ends of the spectrum remote from the transition frequencies of lines of in vivo NMR spectra.
A plot of dispersion vs. absorption mode (DISPA), as described in Sotak et al. Journal of Magnetic Resonance 1984, 57, 453, Herring, F. G.; Phillips, P. S. Journal of Magnetic Resonance 1984, 59, 489, and Craig, E. C.; Marshall, A. G. Journal of Magnetic Resonance 1988, 76, 458, initially used to analyze lineshapes, allows one to determine the phase of an isolated Lorentzian line. If a Lorentzian line is misphased, the DISPA plot will show a circle with diameter equal to the absorption mode peak height rotated around the origin by the number of degrees equal to the phase misadjustment.
Two-dimensional NMR is a well known technique, and is described in Section 8.5 of Alan G. Marshall and Francis R. Verdun, Fourier Transforms in NMR, Optical, and Mass Spectrometry, A User's Handbook, Elsevier, 1990.
For the automated phase correction of symmetrical two-dimensional (2D) NMR spectra with absorptive in-phase diagonal (e.g., NOESY, HOHAHA, z-COSY), Cieslar et al., in Journal of Magnetic Resonance 1988, 79, 154, proposed to maximize the sum of diagonal elements and subsequently to minimize the asymmetry of the diagonal peaks. For 2D spectra with dispersive (e.g., COSY) or absorptive-antiphase diagonal (e.g., DQF-COSY, E. COSY), the reverse must be carried out.
As long as the experimental spectrum has no baseline distortions, infinite signal-to-noise ratio (S/N), infinite digital resolution, and only isolated lines with Lorentzian lineshapes, most methods described will produce a fairly well phased spectrum.
In real spectra, baseline distortions can be caused by several mechanisms and can hardly be avoided. Especially critical are instrumental artifacts like pulse breakthrough and probe ringing that distort the first few points in the FID and cause broad baseline distortions after FT. Long signal averaging times for dilute samples, strong and perhaps incompletely suppressed signals such as solvent lines, aliased dispersive tails of strong signals and unresolved broad resonances are further reasons for baseline distortions. None of the autophasing methods described so far is tolerant concerning these distortions. Especially sensitive are all methods based on integrals and van Vaal's method of determining phases from distant signals by fitting the phase function of a model to both ends of the phase function of the spectrum.
A method capable of phasing spectra at low S/N and low resolution will have to use all the information available concerning the signal phase. Methods based on only a few points around the signal maximum such as DISPA, maximum signal height, and ratio of maximum to minimum signal excursion neglect the phase information in the rest of a signal. Brown's method is limited to the phase information of the baseline and neglects the information of points at the signal maximum.
DISPA plots require at least three points with high S/N around the signal maximum and these points should have magnitudes above 60% of the magnitude maximum. This requires a rather well digitized spectrum. A typical carbon spectrum with linewidths of 0.2 Hz and a spectral width of 200 ppm acquired on a 500 MHz instrument would require a digitization of 2.sup.18 .apprxeq.260,000 complex points, which is significantly higher than normally used for 1D spectra. This might be one of the reasons that Brown found this method to be "extremely unreliable for phasing typical .sup.1 H and .sup.13 C spectra". Similar problems can be expected for all methods involving peak heights unless an interpolation with a lineshape model is involved.
All methods are sensitive to signal overlap. In particular, the maximum of the spectral area doesn't correspond to a well-phased spectrum in the case of signal overlap.
Van Vaals, in U.S. Pat. No. 4,876,507, described a method involving the steps of:
(1) determining peak locations (e.g., for DISPA plots) from a modulus or power spectrum, PA1 (2) correcting overlapping peaks by means of the peak parameters determined, and PA1 (3) using a polynomial as a frequency-dependent phase function determined by a least-squares criterion to determine an overall phase function.
The determination of peak locations from a power spectrum has been shown by Herring and Phillips Journal of Magnetic Resonance 1984, 59, 495, and, using a best-fit polynomial as a phase function has been described by Daubenfeld et al. Journal of Magnetic Resonance 1985, 62, 200, and by Craig and Marshall Journal of Magnetic Resonance 1988, 76, 461.
The Van Vaals patented method, as well as the other known methods described, while providing some measure of phase correction, leave much room for improvement in phase correction methods and results.
Much of the equipment used for obtaining the experimental data relies on a uniform magnetic field or linear magnetic field gradient during excitation of the object to be investigated. Any departure from these standards creates distortions which degrade resolution and signal-to-noise ratio of the data and are very difficult to correct after the data have been acquired. With such distortions it may be difficult or impossible to obtain desired information from the signal. To avoid such distortions, in equipment such as nuclear magnetic resonance equipment, it is important that the magnetic field distortions over the sample volume be as small as possible. Improvements in field uniformity can be made by using additional magnetic fields produced by as many as 38 shim coils placed inside the main magnetic field. Such apparatus is shown in U.S. Pat. No. 3,569,823 as applied to superconducting solenoids and in U.S. Pat. No. 3,622,869 as applied to permanent magnets and electromagnets. The electrical current through these coils may be adjusted to cancel the inhomogeneity of the main magnet, a process called shimming. The application of the invention to magnet shimming is described in U.S. Pat. No. 5,218,299. Other shimming methods are described in U.S. Pat. Nos. 4,987,371 and 5,006,804.
A variety of methods have been proposed for the suppression of solvent signals as described by Tsang et al. in Journal of Magnetic Resonance 1990, 88, 210-215. The proposed methods use experimental signal suppression or correct acquired data through digital filtering techniques or subtraction of a dispersive signal tail function.
For the automated analysis of spectral or imaging data, a criterion to distinguish genuine signals from noise related responses is needed.
Posener, D. W. in Journal of Magnetic Resonance 1974, 14, 121 used parameter precisions to derive equations describing the dependence of the signal detection on the signal lineshape and the spectral S/N ratio. Marshall and Verdun in Fourier Transforms in NMR, Optical, and Mass Spectrometry. A User's Handbook; Elsevier: Amsterdam, 1990; Chapter 5 used parameter precisions to determine the precision of determined peak positions, widths, and heights as a function of peak shape using Monte Carlo type simulations. No applications of various ratios of parameter and corresponding error value (e.g., parameter precisions) for the distinction of genuine signals from noise related responses or for the determination of signal detection probabilities have been used.
Previous methods aimed at using computer-assisted or automated techniques to analyze nD spectral or imaging data have consisted mainly of the application of direct techniques such as cluster analysis and peak picking (Glaser, S.; Kalbitzer, H. R. Journal of Magnetic Resonance 1987, 74, 450 and Meier, B. U.; Madi , Z. L.; Ernst, R. R. Journal of Magnetic Resonance 1987, 74, 565, linear prediction or maximum entropy methods (Gesmar, H.; Led, J. J. Journal of Magnetic Resonance 1989, 83, 53), or least-squares parameter optimization (Madi , Z. L.; Ernst, R. R. Journal of Magnetic Resonance 1988, 79, 513) in cases of severe overlap and/or strong coupling.